Particle kinematics in a dilute, 3-dimensional, vibration-fluidized granular medium 
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We report an experimental study of particle kinematics in a 3-dimensional system of inelastic 
spheres fluidized by intense vibration. The motion of particles in the interior of the medium is 
tracked by high speed video imaging, yielding a spatially-resolved measurement of the velocity 
distribution. The distribution is wider than a Gaussian and broadens continuously with increasing 
volume fraction. The deviations from a Gaussian distribution for this boundary-driven system are 
different in sign and larger in magnitude than predictions for homogeneously driven systems. We 
also find correlations between velocity components which grow with increasing volume fraction. 



The distribution of particle velocities is a fundamental 
descriptor of the statistics of a particulate system. In 
thermal equilibrium, this distribution is always a Boltz- 
mann distribution, but that is typically not the case for 
a nonequilibrium steady state. An archetype of such a 
steady state is the inelastic gas, a system of particles that 
interact by dissipative contact forces, with the energy lost 
being compensated by a driving mechanism. Since the 
inelastic gas is both of fundamental interest and closely 
related to technologically important granular media, it 
has been the subject of much recent experimental, simu- 
lational and theoretical activity [![. Most experiments of 
inelastic gases have focused on 2-dimensional (2D) sys- 
tems; in this article, we present the first detailed measure- 
ment of the velocity distribution in a fully 3-dimensional 
(3D) gas of inelastic grains. 

Inelastic steady states may be broadly divided into 
systems driven from the boundaries and those that 
are driven homogeneously in the bulk. Homogeneously 
heated granular gases are represented theoretically by in- 
elastically colliding particles energized by random, spa- 
tially homogeneous, uncorrelated boosts of energy be- 
tween collisions. It has been shown Q that the high- 
velocity tail of the velocity distribution for this model 
is of the form P(c) tx exp(— Ac 3 / 2 ) where c = v/a is 
a velocity component, v normalized by its r.m.s value 
a. The velocity distribution for low velocities was com- 
puted by finding perturbativc solutions to the Boltzmann 
equation in an expansion in Sonine polynomials around 
a Gaussian [3, In both limits, the results depend on 
inelasticity but not on volume fraction, <fi. These results 
have been confirmed by simulation 0, H, @|. The closest 
experimental analogues for this model are 2D monolay- 
ers of particles fluidized by a vibrated baseQ where ve- 
locity distributions depend strongly on the nature of the 
baseQ ; for a sufficiently rough base , the Sonine expan- 
sion describes the distribution satisfactorily. There are 
no 3D realizations of the homogeneously heated granular 
gas with random forcing in the bulk of the medium. 

Much less theoretical effort has been directed toward 
the more naturally prevalent boundary-driven system. 
On the other hand, several experiments have studied 
2D systems driven by vibration l3, 11 1 and electrostatic 



driving [12j. Some experiments [10L ll2|show the entire 
distribution of velocity fluctuations can be described by 
the functional form P{c) (xexp{~Ac^) with (3 m | over 
a broad range of number density and inelasticity. The 
relationship of these experimental results to the predic- 
tion of Q are unclear since simulations of the homoge- 
neously heated gas [3] show that the asymptotic high- 
velocity behavior only sets in extremely deep in the tail, 
too rare to be experimentally detectable. Furthermore, 
the measured distribution differs between experiments, 
possibly because these quasi-2D experiments are sensitive 
to the specifics of the confinement in the thin dimension 
HQ [Hi or to the substrate on which they move [ll| . 



Three-dimensional steady states have been much less 
studied due to the challenges of tracking particles in the 
interior of a system. Two techniques that have been 
employed for 3D vibrated systems are positron emission 
tracking [l6| of tracer particles, and nuclear magnetic 
resonance imaging [13] • These studies focused on spatial 
profiles of temperature and number density, although a 
non-Gaussian velocity distribution was reported by the 
NMR technique 0. Simulations show mQSjJUllil 
that the velocity distribution evolves continuously from 
a nearly-Gaussian distribution to P(c) ~ exp(—Ac^). 
There are no direct predictions for boundary-driven sys- 
tems, but models 
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231 ] which vary q, the ratio of the 
frequency of particle collisions to the frequency of heat- 
ing events, give some intuition on the passage from the 
homogeneously heated (q » 1) to the boundary-driven 
(g3>l) case. 

In this article we describe experiments using high- 
speed video imaging to directly locate and track particles 
in the bulk of a 3D vibration-fluidized steady state. We 
report kinetic temperature profiles, the form of velocity 
distribution, and comparisons with available predictions. 

As shown in FigQJa), a cubical cell with acrylic verti- 
cal walls is sinusoidally vibrated in the vertical direction 
by an electromechanical shaker (Ling Dynamics V456). 
The vibration frequency used in the experiment ranges 
from 50Hz to 80Hz, velocity amplitudes Vb from 2.3 m/s 
to 3.7 m/s, and accelerations from 90g to 190g. The bot- 
tom and top walls are rough glass plates that provide a 
low-inelasticity surface but also randomise the direction 
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FIG. 1: (a) Sketch of the experimental setup. A cubical cell is 
driven vertically by a shaker. The shaded rectangle is a plane 
of thickness Id illuminated by laser sheets from both sides. 
The darker shading shows the area imaged by the camera, 
(b) Part of a video frame, (c) Edge detection performed on 
this frame, (d) Illustration of the rays emanating from the 
edge, along the intensity gradient, for the top right particle 
in (b) . The solid rays converge toward the center and dashed 
rays diverge, (e) Surface plot for these accumulated gradient 
rays, (f) Located particle positions. For the value of the 
threshold applied to the peaks in (e), the fourth particle in 
the frame is not located. 

of collisional momentum transfer. We use delrin spheres 
of diameter d = 1.560 ± 0.003mm and an average normal 
coefficient of restitution e — 0.92 ± 0.21, experimentally 
determined from particle collisions. The error bar reflects 
the dependence of e on impact parameter, relative par- 
ticle velocity, and spin. The side of the cell is 51.2mm 
« 32<i and is illuminated by light sheets produced by 
expanding beams from laser diodes (Thorlabs ML101J8) 
with a cylindrical lens. The thickness of the sheet in 
the y-direction is « Id, and its y-position can be varied, 
allowing us to study particle motions in x-z planes at 
varying depth from the front wall. A Phantom v7 cam- 
era images the x-z plane selected by the light sheet, at 
5,000 frames/second, with a resolution of 640 x 480 pix- 
els. The field of view is a rectangle of width 30d and 
height 20c?, centred on the middle of the cell. 

The difficulty in locating particles stems from two dis- 
tinct eclipsing effects. First, particles in the middle of the 
cell are less likely to be illuminated because the light- 
sheet is obstructed by particles in the light path. The 
delrin spheres are homogeneously illuminated; any laser 
light incident on a sphere is scattered through its vol- 
ume. However, when the light-sheet is at some depth 
from the front plane, illuminated spheres can be partially 
eclipsed by particles in front of them. We focus on ve- 
locity statistics rather than number density distributions 
since the detection probability of a particle depends on 
these two effects. A video frame with examples of par- 
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FIG. 2: Distribution of the kinetic temperature at planes of 
various j/(top row) and z(lower row). The dashed line shows 
the midplane of the cell. The driving frequency is 60Hz, T — 
HOg and volume fraction 5.1%. 

tially eclipsed particles is shown in FigQ}). To locate 
particles, we exploit the fact that the convex illuminated 
edges of the images are circular arcs of known radius. 
We find illuminated edges (FigJTh) and draw rays from 
all edge points in the direction of the local intensity gra- 
dient. For points on convex illuminated edges, these rays 
converge to the centre, while concave edges produce rays 
that diverge (FigJTJi). Local maxima for accumulations 
of rays (FigQJ) above a cut-off value are candidates for 
particle centres; objects with an insufficient length of il- 
luminated perimeter are rejected. Particle centres are 
determined with subpixel resolution by minimising the 
squared distance to the gradient rays. Eclipsing produces 
no systematic bias in locating particle centres. 

The kinetic temperature T — \{{v x ) + (v%)) i s mea- 
sured at several depths, y, between the front plane and 
a little beyond the middle of the cell. As shown in Fig[2] 
(top row), in any given x-z plane, T is higher near the top 
and bottom wall, and lower in the middle. Comparing 
x-z planes at different depths y, we see that T goes up 
in the middle of the cell, revealing the dissipative effect 
of the front vertical wall. The data are then reorganized 
to show T(x,y) at a few heights, z (Fig|2] bottom row). 
Within each plane, T is lower close to the vertical walls, 
again due to the dissipation at the walls. The length 
scale over which the wall dissipation manifests itself is 
comparable to the mean free path, and grows at lower <f>. 

We now turn to the entire distribution of the velocity 
fluctuations. The velocity fluctuations are anisotropic; 
we concentrate on v x , the horizontal component, perpen- 
dicular to the driving direction. The distribution P(c x ) of 
c x , the normalized horizontal velocity c x — v x /(2T x ) 1 ^ 2 
where T x = (v 2 ), shows weak dependence on position 
very close to the walls of the cell (as also seen in sim- 
ulations 0, H3|). Therefore, we report data at y = 8d 
which is far from the front wall, and yet not so deep in 
the cell where the velocity statistics are greatly dimin- 
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FIG. 3: P(c x ) for (f> — 5.1% for various driving frequencies 
and amplitudes, sampled at y = 8d. One data set taken with 
smooth top and bottom plates, and another set taken at y = 
3d are included for comparison. The inset shows the average 
kinetic temperature of the sampled region as a function of Vq, 
where Vo is the amplitude of the driving velocity. 





FIG. 4: P(c x ), different volume fractions ranging from 1.8% 
to 7.2%. The smooth lines connecting symbols are a guide to 
the eye. The inset shows ln(~ln(c x )) vs ln(c x ). 

ished by eclipsing effects. The dependence of P(c x ) on 
different driving frequencies and amplitude is shown in 
FigEl As the overall temperature is changed by a fac- 
tor of 3, we observe no systematic changes in P(c x ). In 
experiments on 2D monolayers @, 0] it was noted that 
velocity distribution depended on the smoothness of the 
driving surface, as might be expected when interparti- 
cle collisions and heating events occur with comparable 
frequency. To test whether the influence of the bound- 
ary persists into the interior, we replace the rough glass 
plates by smooth delrin plates and do not observe any 
change in P(c x ). This suggests that the observed statis- 
tics are a consequence of inter-particle collisions, and are 
insensitive to the details of the driving surface. 

Having established that distribution of horizontal ve- 
locities is insensitive to T, to the driving surface, and to 
location within the cell, we now discuss the dependence of 
P(c x ) on the number density of particles. We have stud- 
ied six volume fractions cf> ranging from 1.8 to 7.2%. The 
lower limit on 4> is chosen so that the mean free path is 
still much smaller than the cell dimensions, and the upper 



FIG. 5: (a)P(cx) at volume fractions <f> = 1.8,3.9 and 6.2% 
fitted to second-order Sonine polynomials with 02 treated as 
a free parameter. (b)P(c x ) at volume fraction <f> — 6.2% with 
best fit second- (red dash line, 02 = 0.440) and third-order 
(red dotted line, 02 = 0.455, a 3 = -0.0287) Sonine forms. The 
green dashed and dotted lines show the 2nd and 3rd order 
Sonine forms with the coefficients 02 and 03 computed from 
the measured cumulants. (c)The deviation function D(c x ) at 
small c x for the above volume fractions. (d)The best fit value 
of ai (a) compared with 0,2 (□) computed from (c 4 ). 



limit is constrained by poor detection statistics at high 
volume fractions. Deviations from a Gaussian are appar- 
ent even at the lowest </>. With increasing (j>, the tails of 
P(c x ) get broader. Thus the velocity distribution varies 
continuously with density unlike in some 2D experiments 
[13, E3] where P(c x ) is unchanged over a broad range of 
<p. This is also unlike predictions for the homogeneously 
driven or cooled state [2| where P(c) is independent of </>. 

The high-velocity tail of P{c x ) cannot be described 
by the form exp(-cP): as shown in the inset of FigJJJ 
a plot of ln(—lnP(c x )) against ln(c x ), shows curvature, 
whereas in the equivalent 2D experiment (loj we observed 
a straight line with a slope (3 — 1.55 ± 0.1. The statistics 
in the experiment only capture the tail up to 4 decades 
below the peak, and leave open the possibility that this 
could be the asymptotic form of the distribution at large 
c. However, for any realistic description of grain dynam- 
ics, even rarer fluctuations are probably irrelevant. Ear- 
lier simulations [H, [H, [H [U [1| have found density- 
dependent velocity distributions, however, this is the first 
experimental study in 3D to observe this effect. 

In the absence of predictions for the a boundary-driven 
system, we compare P(c) to predictions for a homoge- 
neously heated inelastic gas, where at small velocities, 
the deviations D{c) from the gaussian distribution, $(c), 
have been perturbatively calculated 0, Q as an expansion 
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FIG. 6: The correlation of v x and v z (•), and the correlation 
implicit in the Sonine formula(A): 1 + d2- The inset shows 
the anisotropy, (v\ — v^)/(v x + v^). 



in Sonine polynomials S p (c 2 ) : 
P(c) = $(c) [1 + D(c)] = $(c) 



1 + ^ apS'p 



(1) 



i(i and 



The first two polynomials are S\{c 2 ) = —c 
S 2 (c 2 ) = ±c 4 -§(<? + 2)c 2 + \d{d + 2) where d is the 
dimensionality. The coefficients a p are given in terms 
of the moments of the distribution, P(c): cii = 0, a 2 — 
4(c 4 ) /(d(d + 2)) - 1. Predictions for the dependence 
of a 2 on the restitution coefficient for the homogeneously 
heated and cooling states have been validated by both 
DSMC and event-driven simulations [l, 0, In FigOwe 
test whether the Sonine expansion is a good description 
of P(c x ) in our boundary driven system by fitting the 
data for varying volume fraction to a second-order Sonine 
expansion with a 2 as a free parameter. 

FiglH a )"(c) show that the second-order Sonine correc- 
tion works best at the lowest volume fraction and the 
range over which D(c), the deviation from a gaussian, 
is well-fit diminishes at higher volume fractions. As in 
the quasi-2D experiment of which models a homoge- 
neously heated gas; the quality of the fit is reasonable for 
c < 2. A third-order Sonine term does not improve the fit, 
as shown in FigOJb). 

Apart from the dependence on <j) of the fit param- 
eter a 2l we also note that a 2 is opposite in sign, and 
much larger in magnitude than that found in the homo- 
geneously heated or cooled states for the same nominal 
restitution coefficient. Furthermore, as shown in FigJS^d) 
the best-fit value of a 2 disagrees with the value directly 
calculated from the 4th cumulant, (c 4 ), raising the possi- 
bility that the gaussian reference state may not be appro- 
priate for a boundary-driven system. To our knowledge, 
the only predictions for P(c) in a boundary-driven sys- 
tem were made in 24j , where the 4th cumulant is treated 
as an independent hydrodynamic field. They find (c 4 ) 
shows a density dependence qualitatively like ours, but 
with much lower magnitude than we measure. 

A general consideration regarding P(c) is that if a ve- 
locity distribution is isotropic and if the velocity compo- 



nents are uncorrelated, it must be a gaussian distribution. 
A non-gaussian distribution implies that one or both of 
these assumptions is invalid. For low volume fractions, 
the velocity fluctuations are anisotropic, but they become 
more isotropic at higher densities. This trend is shown 
in the inset to FigUJ where the anisotropy is quantified 
by 2 (v 2 — v 2 ) j (v 2 + v 2 ) . However, since the velocity 
distribution does not tend towards a Gaussian at large 
</>, the velocity components must be correlated as shown 
in Fig|6] Indeed, the correlation between velocity com- 
ponents (v 2 v 2 ) j (y 2 ) (v 2 ) grows with volume fraction </>. 

Our measurements of the 3D particle kinematics in the 
interior of a vibration-fluidized granular medium thus re- 
veal a non-Gaussian velocity distribution that is insensi- 
tive to conditions at the driving surface. The shape of 
the distribution evolves continuously with volume frac- 
tion; the functional form differs markedly from the ho- 
mogeneously heated state, thus emphasizing the need for 
theoretical development for boundary-driven systems. 
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